% bci2ecog
% bciGetDataSample
% ecogShowVariance %where is energy
% ecogPCS %show PCs, allow for feature selection 
% optimize spectrogram estimation
% notch filtering (AR filter ??)

%% get data and set parameters
load EMG256Ch3
load EMG256Ch3_Epoch
load EMG256Ch3_Rawdata
ecog.data=single(trainDat);
ecog.sampDur=1000/bci.TDT.Fs(3);
ecog.selectedChannels=1:2;
ecog.nSamp=size(ecog.data,2);
ecog.nBaselineSamp=0;
ecog.timebase=single(0:ecog.sampDur/1000:ecog.sampDur*(ecog.nSamp-1)/1000);
clear trainDat trainEpoch
%Downsample timeseries
ecogDS=ecogDownsampleTS(ecog,200);
dsFactor=ecogDS.sampDur/ecog.sampDur;
% OPTIONAL: remove common average reverence
%ecog=ecogRemoveCommonAverageReference(ecog);
% Make Segments
windowOnsetIdx=Dat.trainOnset; 
windowLength=bci.init.nWinLength; % that's in samples
windowLabel=Dat.trainLabel;
ecogSeg=ecog;
[ecogSeg]=ecogSegmentTS(ecogSeg,windowOnsetIdx,0,windowLength);
ecogSeg.data=single(ecogSeg.data);
ecogSeg1=ecogSeg;
ecogSeg1.data=single(ecogSeg.data(:,:,find(windowLabel==1)));
ecogSeg2=ecogSeg;
ecogSeg2.data=ecogSeg.data(:,:,find(windowLabel==2));
ecogSeg3=ecogSeg;
ecogSeg3.data=ecogSeg.data(:,:,find(windowLabel==3));
%% Check 1: visual inspection 
ecogTSGUI(ecog)
%% Check 2: Periodograms show bad frequency ranges 
params.Fs=1000/ecog.sampDur;
params.fpass=[1 100]
params.tapers=[4 7];
ecog=ecogMkPeriodogramMultitaper(ecog,1,params)
ecogP=ecogRaw2Ecog(log(1e20*ecog.periodogram.periodogram)',0,1000*diff(ecog.periodogram.centerFrequency(1:2)),[]);
ecogTSGUI(ecogP)
ecog=ecogFilterTemporal(ecog,[66 55],[3 3]);
%% Check 3: standard deviation IMPORTANT FOR ARTIFACT THRESHOLS
figure;plot(std(ecog.data'))
%% Check 4: spectrogram 
params.Fs=1000/ecog.sampDur;
params.fpass=[1 100]
params.tapers=[10 9];
params.trialave=0; 
ecogDS=ecogMkSpectrograms(ecogDS,1,params) 
for k=1:length(ecogDS.selectedChannels)
    figure;imagesc(squeeze(ecogDS.spectrogram.spectrogram(:,:,ecogDS.selectedChannels(k))))
    title(['Channel ' num2str(ecogDS.selectedChannels(k))])
end
max(max(ecogDS.spectrogram.spectrogram(:,:,1)))/10
set(gca,'clim',[-1.7337e-018 1.7337e-018])
%% Check for artifacts

%% Check 5: periodogram on time series segments
[ecogSeg.data,ecog.segmenOffset]=ecogSegmentTS(ecog.dat,windowOnsetIdx,0,windowLength);
ecogSeg1=ecogSeg;
ecogSeg1.data=ecogSeg.data(:,:,find(windowLabel==1));
ecogSeg2=ecogSeg;
ecogSeg2.data=ecogSeg.data(:,:,find(windowLabel==2));
ecogSeg3=ecogSeg;
ecogSeg3.data=ecogSeg.data(:,:,find(windowLabel==3));
ecogSeg=ecogMkPeriodogramMultitaper(ecogSeg,1,params);
ecogSeg1=ecogMkPeriodogramMultitaper(ecogSeg1,1,params);
ecogSeg2=ecogMkPeriodogramMultitaper(ecogSeg2,1,params);
ecogSeg3=ecogMkPeriodogramMultitaper(ecogSeg3,1,params);
% periodogram to data field
ecogSegP1=ecogRaw2Ecog(log(1e20*permute(ecogSeg1.periodogram.periodogram,[2 1 3])),0,1000*diff(ecogSeg1.periodogram.centerFrequency(1:2)),[]);
ecogSegP2=ecogRaw2Ecog(log(1e20*permute(ecogSeg2.periodogram.periodogram,[2 1 3])),0,1000*diff(ecogSeg2.periodogram.centerFrequency(1:2)),[]);
ecogSegP3=ecogRaw2Ecog(log(1e20*permute(ecogSeg3.periodogram.periodogram,[2 1 3])),0,1000*diff(ecogSeg3.periodogram.centerFrequency(1:2)),[]);
% get t-values
pT1Vs2=ecogTValues(ecogSegP1,ecogSegP2);
pT3Vs2=ecogTValues(ecogSegP3,ecogSegP2);
pT1Vs3=ecogTValues(ecogSegP1,ecogSegP3);
figure; plot(ecogSegP1.timebase/1000,pT1Vs2);xlabel('Frequency [Hz]');title('1 vs 2');ylabel('t-value');legend(num2str([1:size(ecogSegP1.data,1)]'))
figure; plot(ecogSegP1.timebase/1000,pT3Vs2);xlabel('Frequency [Hz]');title('3 vs 2');ylabel('t-value');legend(num2str([1:size(ecogSegP1.data,1)]'))
figure; plot(ecogSegP1.timebase/1000,pT1Vs3);xlabel('Frequency [Hz]');title('1 vs 3');ylabel('t-value');legend(num2str([1:size(ecogSegP1.data,1)]'))


% check 5: trial average spectrogram
%% check classifier
bciLoad('EMG256Ch',1,1);
bciRefreshParam;
bciRecreate(1);
bciCheckClassifier(1);
